Supersymmetry solution for finitely extensible dumbbell model 
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Exact relaxation times and eigenfunctions 
for a simple mechanical model of polymer dy- 
namics are obtained using supersymmetry meth- 
ods of quantum mechanics. The model includes 
the finite extensibility of the molecule and does 
not make use of the self-consistently averaging 
approximation. The finite extensibility reduces 
the relaxation times when compared to a linear 
force. The linear viscoelastic behaviour is ob- 
tained in the form of the "generalized Maxwell 
model". Using these results, a numerical inte- 
gration scheme is proposed in the presence of a 
given flow kinematics. 

83.10.Nn Polymer dynamics 
05.10.Gg Stochastic analysis methods 
03.65.Fd Algebraic methods 

I. INTRODUCTION 



Simple mechanical models are helpful for 
understanding the complex dynamical be- 
haviour of polymer molecules. The so-called 
elastic dumbbell constitutes the simplest 
model that captures the effects of stretch- 
ing and orientation of the polymer In 
this model the polymer is represented by two 
beads which are connected by a spring. The 
classical Hookean dumbbell model |l[] is char- 
acterised by a linear spring force. Its relax- 
ation times decrease like the inverse of the 
mode number and the relaxation modulus is 
of a single exponential form. As noted by 
many authors, the finite extensibility of the 
molecules has to be included in order to avoid 
unphysical behaviour such as infinite viscosi- 
ties predicted by the Hookean model. It is 
even argued "that taking into account the 
nonlinear stretching of the polymer molecules 
is the most important correction that should 
be made to the Hookean model in order to 
describe real systems" (Ref. [Q, Chap. 13, p. 
81). Molecular arguments suggest the finitely 
extensible spring force law to be given by the 
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inverse Langevin function Very often, ra- 
tional H or Pade approximations [|| to the 
inverse Langevin function are used since they 
are analytically more tractable. All these ap- 
proximations show the desired physical char- 
acteristics: a linear regime for small exten- 
sions of the spring and a divergence at a fi- 
nite extension. These forces therefore pro- 
hibit the molecule to stretch beyond a certain 
value which is related to the total length of 
the chain. 

In general, the finite extensibility of the 
molecules prohibits analytical expressions for 
the spectrum of relaxation times and the re- 
laxation modulus. The shear relaxation mod- 
ulus for the Warner spring force (so-called 
FENE model, jjJD) nas been determined nu- 
merically Q. The results are very close to 
those obtained for the Peterlin approxima- 
tion to the Warner force (FENE-P model, 
P]||). The effect of the finite extensibility 
on the spectrum of relaxation times remains 
less studied. One might speculate that the 
finite extensibility decreases the slowest re- 
laxation time as described by the Peterlin 
approximation. However, knowledge of the 
whole spectrum of relaxation times would be 
very desirable for understanding the dynam- 
ics of dumbbell models. It would allow, for 
example, to identify the number of relevant 
slow modes that contribute to the polymer 
contribution to the stress. Work in this di- 
rection was performed in M where the univer- 
sal constitutive equation of dilute polymeric 
solutions was found. There, it was assumed 
that in the limit of low Deborah and Weis- 
senberg number the modes corresponding to 
the two lowest eigenvalues dominate. 

Here, we compute the relaxation times, 
eigenfunctions and the relaxation modulus 
of a finitely extensible dumbbell model ex- 
actly. The model is a one-dimensional ver- 
sion of the finitely extensible dumbbell mod- 
els that uses trigonometric functions to ap- 
proximate the inverse Langevin force. It is 
demonstrated that this force approximates 
the inverse Langevin function as good or 
even better than the Warner force frequently 
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used in Brownian dynamics simulations. The 
present study therefore clarifies the effect of 
the finitely extensible spring force on the 
spectrum of relaxation times, eigenfunctions 
and the relaxation modulus. 



II. MODEL DESCRIPTION 

Let f(q;t) denote the distribution function 
of the connector vector q of the dumbbell at 
time t. The time evolution of / is given by 
the well-known Fokker-Planck equation [jD 

d t f = DV{{VU]+V}f = L FP f (1) 

where we have assumed that the spring force 
F can be derived from the dimensionless po- 
tential U via F(q) = k B TVU(q). The dif- 
fusion constant D = 2k B T/C, is assumed to 
be independent of q which means we neglect 
hydrodynamic interactions. Note that the 
equilibrium distribution / eq is given by / oq = 
N oq e~ u , N oq being a normalization constant. 
As is well-known the Fokker-Planck equation 
(Q) can be transformed into an equivalent 
imaginary time Schrodinger equation ||^| . If 
ip is a solution to 

- d t ip = DHyj (2) 

with the 
operator H = -{l/D)e u / 2 L FP e~ u / 2 , then 

1/2 

the function f(q;t) = ip{q\t)f C q (q) solves 
the original Fokker-Planck equation ([!]). In 
particular, solving the eigenvalue problem 
Hip n — e n ip n leads immediately to the so- 
lution of the eigenvalue problem of the origi- 
nal Fokker-Planck operator Lfp/Vi = —K^fn 
with A^ 1 = De n and /„ = ip n fcq 2 , where fac- 
tors D are introduced for later convenience. 
Inserting (Q) for the Fokker-Planck operator 
H yields the form of the Schrodinger operator 

H = -V 2 + V{q) 1 (3) 

where the Schrodinger potential V is given by 
V = $ 2 - V$, with $ = VC//2. Therefore, 
each solvable potential for the Schrodinger 
equation serves as a solvable model for the 
Fokker-Planck equation. In the following, we 
benefit from the fact that simple Schrodinger 
potentials V can have quite complicated 
counterparts U for the Fokker-Planck equa- 
tion. Note also that operator H is Hermitian 
which is in general not the case for Lfp- 



The simplest exactly solvable potential for 
the Schrodinger equation related to finite ex- 
tensibility is the three dimensional spheri- 
cal symmetric infinite well potential. Up to 
a constant, the corresponding potential U 
of the Fokker-Planck equation is given by 
U = — 2kn/> cq , where ip eq is the ground state 
of the Schrodinger problem. In the present 
case, the ground state is known to be tp eq ~ 
sm(irq/qo)/q with q = \q\ leading to the force 
law 

F(q) = 2k B T 1 - { * q/qo) q COt{nq/q ° ) l. (4) 

Since the force (0) is one-parametric, the 
"spring constant" is no adjustable parame- 
ter but determined by the maximum exten- 
sion q . If we define the spring constant h in 
the linear regime of (^), the dimensionless fi- 
nite extensibility parameter b — hq^/(k B T) is 
fixed: b = 63 = 27r 2 /3. Parameter b is known 
to denote roughly the number of monomer 
units represented by the dumbbell and should 
therefore be a large number. Thus, (Q) is of 
limited use. However, if b = 63 is acceptable 
the force (Q) compares very well to the inverse 
Langevin force and Cohen's Pade approxima- 
tion, as shown in the inset of Fig. |l|. In this 
case, exact eigenfunctions of the Schrodinger 
operator are known to be the spherical Bessel 
functions and the exact eigenvalues are given 
by the zeros of this functions. 

Let us now consider the one-dimensional 
case. One-dimensional nonlinear dumb- 
bell models serve as toy models for vari- 
ous approximations that have been proposed 
to obtain closed constitutive equations for 
polymer solutions jl0|]. In addition, the 
one-dimensional models are believed to de- 
scribe the elongational behaviour of the poly- 
mer molecule reasonably The one- 
dimensional version of the operators L-pp and 
H is obtained by replacing q with q and V 
with d/dq. 

To account for the finite extensibility we 
propose the following two parameter family 
of force laws 

F(q) = hqi tan(q/qi), for - q < q < q Q 

(5) 

where q\ = 2q$ /it is determined by the maxi- 
mum extension qo and h denotes the "spring 
constant" in the sense F(q) — > hq for q — > 0. 
In Fig. [l] we plot the force F/h as a func- 
tion of the reduced extension q/qo- Fig. 
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|l| also shows the Warner force || F(q) = 
hq/(l — g 2 /g 2 ), the inverse Langevin force 
(h/3)L- 1 (q/q ), with L(x) = coth(x) - 1/x 
and Cohen's Pade approximation F(q) = 
(hq/3)(3 - <7 2 /<7o)/(l - Q 2 /l 2 ) §■ It is eas- 
ily verified that these approximations give the 
correct limiting behaviour not only for small 
but also near the maximum extension. In ad- 
dition, as can be seen from Fig. [I], the force 
serves as an even better approximation to 
the inverse Langevin force as does the Warner 
force. Therefore we accept as a reason- 
able force law to substitute for the inverse 
Langevin force in one dimension. 

The potential U from which the force (||) 
can be derived is U{q) = — b\ lncos(g/gi), 
where we introduced the dimensionless pa- 
rameter bi = hq\/(k B T) = (2/tt) 2 &. 
The equilibrium distribution function of 
the Fokker-Planck equation is f e q(q) = 
N cq cos bl (q/ qi ), with jV^ 1 = q X 2^B[{bx + 
l)/2, (&i + l)/2], where B[x,y] is the Beta 
function. The equilibrium distribution is very 
close to the one corresponding to the Warner 
force f cq (q) = N cq (l~q 2 /q 2 ) b / 2 , with iV^ 1 = 
gbfl[l/2,(6 + 2)/2]. 

III. SUPERS YMMETRY SOLUTION 

To obtain exact eigenvalues and eigenfunc- 
tions of the Fokker-Planck equation (Q) we 
exploit supersymmetry methods for the cor- 
responding Schrodinger problem. Note, that 
the model force (H) belongs to the class of 
so-called Poeschl- Teller potentials which may 
also be solved exactly by Schrodinger's fac- 
torization method p2fl . Consider the orig- 
inal potential U together with the inverted 
potential — U and denote by /_ the solution 
for the potential U and / + for —U. As done 
before, defining f± — ^N cq e ±u ^ 2 ip± we ar- 
rive at — (l/D)dtip± = H±ip± with opera- 
tors H± = V 2 + V±(q). The "partner po- 
tentials" are given by V± (q) = $ 2 (q) ± $' (q) 
with $ = U'/2, where $' and U' denote the 
derivative of $ resp. U with respect to q. 
Schrodinger operators H± are very well stud- 
ied in the literature since they occur in the 
so-called "Witten model" which is the 

simplest model that shows all typical features 
of supersymmctric quantum mechanics. 

Rewriting the force (|1) as F = 2k^T<^ 
we obtain <f> a (q) = (xtaxi(q/q{), for —qo < 
q < qo- We denote explicitly the depen- 
dence on the parameter a = b\/(2q{) which 



is proportional to the spring constant h. To 
proceed further we take advantage of the 
concept of "shape invariance" introduced by 
Gendenshtein The partner potentials 

V±{a,q) = ± $' a (q) are called shape 

invariant if they are related by V+ (ak , q) = 
V- (afc+i , q) + R{au+i ) where the new param- 
eter otu+x is a function of a^. In our case, 
$ Q = atan(g/gi), shape invariance is easily 
verified for otk+i = Oik + ao = a and 

R(ak+i) — <^k+i ~ a k- F° r shape invariant 
potentials it is possible to define a series of 
operators that are isospectral except for the 
lowest eigenvalue which is J2k=i R( a k)- Go- 
ing back to the original Schrodinger opera- 
tor with lowest eigenvalue zero the complete 
spectrum is found to be e„ = X)fe=i R( a k) 
0,0. Inserting the special form of R and 
multiplying by D, the exact (inverse) eigen- 
values of the Fokker-Planck operator (Q) are 
given by 

A„ = 2A H [n + n 2 /&i]~\ n = l,2,... (6) 

where Ah = C/(^h) is the time constant of the 
Hookcan dumbbell. Note, that the eigenvalue 
zero corresponds to the equilibrium distribu- 
tion. Corresponding eigenfunctions of opera- 
tor H are 

Mq) = C n a f a\ . . . 4_i cos n+b ^ 2 (q/ qi ) (7) 

where a k = —d/dq + <& ak are generalized 
creation operators and C n denote normaliza- 
tion constants |l2). Remember that eigen- 
functions of the original Fokker-Planck op- 
erator are obtained by /„ = V-Vi/cq 2 - The 
first eigenfunctions read: /o = f cqi fi = 
Ni sm{q/ qi )f cq , h = N 2 [(h + 1) sm 2 (q/ qi ) - 
cos 2 (q/qi)]f cq with Nf = h + 2 and iVf = 
(6i+4)/[2(6i + 1)]. Note, that the eigenfunc- 
tions are orthonormal, (/„, f m ) — S nm , in the 
scalar product 

/qo 
f- 1 (q)g(q)h(q)dq. (8) 
-90 

The linear viscoelastic behaviour can be ob- 
tained from linear response theory. For a 
given flow kinematics the perturbation of the 
Fokker-Planck operator is given by L cxt (i) = 
—d q n(t)q where ft(t) denotes the velocity gra- 
dient. 

The polymer contribution to the stress, 
r/(fc B T) = (a J) -l,a = f eq qU', is found 
to be 
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r/(k B T) = G(t- t')K{t')dt', (9) 

J —OG 

where the relaxation modulus G(t) is given by 
the equilibrium correlation function G(t) = 
(a,e LFpt <r) — 1. Inserting the completeness 
relation of the eigenfunctions, S(q — q') = 

EnMQ)Ml') = /cqV)E„/n(9)/n(9'), 

we obtain 



G(t) 



OO 

E 

n=l 
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which is of the form of the "generalized 
Maxwell model" . Note, that only even terms 
are included in ( |To|) since they correspond to 
symmetric eigenfunctions whereas odd terms, 
corresponding to antisymmetric eigenfunc- 
tions, vanish by symmetry. Analytical ex- 
pressions for G(t) resulting from dumbbell 
models are known for linear springs (Hookean 
dumbbell) where G[t) = e~*/ AH . For the 
Warner spring force no such expression is 
available. The Peterlin approximation (!]] 
leads also to a single-exponential form with 
relaxation time Ah [1 + 1 /b] ~ 2 . The relaxation 
modulus ([l0]) in our model shows a spectrum 
of relaxation times A„, the lowest one be- 
ing A 2 = A H [1 + 2/fri]- 1 . The finite exten- 
sibility of the dumbbell therefore reduces the 
longest relaxation time when compared to the 
Hookean dumbbell. Only in the limit b — » oo 
the result for the Hookean dumbbell is recov- 
ered. As can be seen from Fig. g G(t) is dom- 
inated by the lowest relaxation time since the 
relative weight of the higher modes decreases 
rapidly. For b — 20, the relative weight of 
the first two modes is 99%, for b — 50 the 
first mode carries 99% and for b — > oo all 
weights except for the first vanish, reflecting 
the single-exponential form of the Hookean 
dumbbell. 

The relaxation modulus determines the 
zero-elongational viscosity f] — J °° G{t)dt. 
As shown in , 770 can be expressed in terms 
of second moments of the equilibrium distri- 
bution function for arbitrary dumbbell mod- 
els. In our model, the resulting integral can- 
not be done analytically. 

Knowledge of the exact relaxation times 
and eigenfunctions (||) and (Q) can further be 
used to integrate the kinetic equation (|l|) in 
the presence of a given flow field. Expanding 
the distribution function into eigenfunctions, 
/fe*) = E„ c « (t)fn(q), leads to 
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with A n k = (f n , d q qfk) a constant matrix in- 
dependent of the flow. Coefficients c n (t) are 
defined by c n (t) = (/„,/(«; t)>. Eqs. ^ 
are an equivalent formulation of the kinetic 
equation (|l|). Conservation of total proba- 
bility is guaranteed by cq = 0, which is eas- 
ily verified from (p"l|). The stress r is deter- 
mined by the coefficients c„ via T/(k-g,T) = 
Sr=i( cr '/2r l }c2n- Because of symmetry, only 
even terms contribute to r. Moreover, Eqs. 
( p"T| ) for even coefficients decouple from those 
for odd coefficients since A n k is nonzero only 
if n and fc are both even or both odd. There- 
fore, if we are interested only in the polymer 
contribution to the stress, only even coeffi- 
cients have to be included in ([fl]) . Given A n k , 
Eqs. ( |Tl| ) can be integrated for a finite set of 
eigenfunctions quite efficiently using standard 
integrators for ordinary differential equations. 
How many eigenfunctions have to be included 
in order to obtain a desired accuracy is a del- 
icate question. A simple estimation is the 
following. To a first approximation, coeffi- 
cients c„ are still of exponential form but 
with a different time constant A^ 1 + nA nn . 
Modes that decrease quadratically with num- 
ber n will then be negligible. This leads to 
a minimal number of modes n c that satisfies 
n^/bi = n c + 2\n\K,(t)A ncnc |. If this number 
is not too large, method (jll]) is much more 
efficient than Brownian dynamics simulation. 



IV. CONCLUSIONS 

We have obtained the spectrum of relax- 
ation times, the eigenfunctions and the relax- 
ation modulus for a one-dimensional finitely 
extensible dumbbell model exactly. It is ob- 
served that the relaxation times are reduced 
when compared to the Hookean dumbbell 
model. Moreover, the decrease of the re- 
laxation times with mode number n shows 
a crossover from vT 1 to nr" 1 at the value of 
the finite extensibility parameter b. The re- 
laxation modulus is obtained in the form of 
the "generalized Maxwell model" . For rather 
large values of the finite extensibility param- 
eter 6, which are commonly believed to be 
physically meaningful (see, e.g. (|)), the lin- 
ear viscoelastic behaviour is well described by 
a single relaxation time. This is because the 
relative weight of the higher modes is decreas- 
ing rapidly with increasing b. Therefore, it 
is not surprising that the linear viscoelastic 
regime is well described by the Peterlin ap- 
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proximation as found in lq|. To calculate the 
non-linear dynamical response to a given flow 
kinematics we propose a numerical integra- 
tion scheme using the expansion of the distri- 
bution function into the exact eigenfunctions 
in the absence of flow. 

All these results are limited to the one- 
dimensional case. In three dimensions the 
class of solvable potentials is much more re- 
stricted. In particular we did not find any 
relevant two parameter potential that can be 
solved exactly. However, there is some evi- 
dence, that "a simple one-dimensional version 
of the FENE theory captures qualitatively 
and quantitatively the clongational behaviour 
of the actual three-dimensional theory" (Ref. 
[fiof , Introduction). 
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FIG. 1. Nonlinear, finitely extensible dumb- 
bell forces are shown as a function of the re- 
duced extension. Long dashed: Warner force, 
dotted: the force (|E]) proposed here, dashed: Co- 
hen's Pade approximation, solid line: the inverse 
Langevin function and dot-dashed: the Hookean 
spring. As easily seen from the figure the proposi- 
ton (Eh approximates the inverse Langevin func- 
tion even better than does the most frequently 
used Warner approximation. The inset shows 
these forces for the special case b = 63 in three di- 
mensions. The same symbols are used as before, 
but the force (@) is replaced by (0). 
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FIG. 2. The first ten relaxation times X n , 
n even, are displayed as full circles (^j). 
The histogram shows their relative weights 
\{a, f n )\ 2 /G(0). The dimensionless parameter b 
was chosen to be b — 20. With increasing b the 
relaxation times approach 1/n while A2 accumu- 
lates the total weight. 
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